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We present a detailed numerical study of effective interactions between micron-sized silica spheres, 
induced by highly charged zirconia nanoparticles. It is demonstrated that the effective interactions 
are consistent with a recently discovered mechanism for colloidal stabilization. In accordance with 
the experimental observations, small nanoparticle concentrations induce an effective repulsion that 
counteracts the intrinsic van der Waals attraction between the colloids and thus stabilizes the 
suspension. At higher nanoparticle concentrations an attractive potential is recovered, resulting 
in reentrant gelation. Monte Carlo simulations of this highly size-asymmetric mixture are made 
■ possible by means of a geometric cluster Monte Carlo algorithm. A comparison is made to results 

' obtained from the Ornstein-Zernike equations with the hypernetted-chain closure. 
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PACS numbers: 61.20.Ja, 64.75. +g, 82.70.Dd 



I. INTRODUCTION 



Colloidal suspensions are used in a wide range of applications, such as coatings and drug carriers, and as a precursor 
for various advanced materials, including colloidal crystals. In all these systems, the stability of the suspension 
plays a central role. The van der Waals attractions between the colloids — if not counteracted by some appropriate 
mechanism — lead to their aggregation, typically resulting in the formation of a gel. Conventional methods to mitigate 
q the effect of the van der Waals forces rely on charge stabilization or steric stabilization [lj . In charge stabilization, an 
l—— electrostatic repulsion between the colloids is induced, e.g., by an appropriate choice of the pH. Steric stabilization 
frequently involves grafting short polymer chains onto the colloidal surface, which prevents a close approach of the 
colloidal surfaces. 

■<^j- ■ Quite recently, an alternative strategy for colloidal stabilization has been discovered Q . These experiments involve 
an aqueous suspension of silica spheres with a diameter <7 m j cro of approximately 0.6/xm, to which charged zirconia 
nanoparticles of much smaller diameter (<7 nano — 6nm) are added. In the absence of nanoparticles, as well as in the 
presence of very small concentrations of nanoparticles, the microspheres exhibit a tendency to aggregate owing to 
their van der Waals attraction. Increasing the nanoparticle concentration prevents this aggregation and thus stabilizes 
the suspension. At even higher zirconia concentrations, the aggregation behavior reappears, leading to a "window 
of stability" in the nanoparticle concentration. In order to ensure that the observed stability does not result from 
direct electrostatic repulsion between the microspheres, the suspension is kept at a very low pH, near the isoelectric 
point of silica. Furthermore, it has been verified in Ref. |2| that there is no strong adsorption of nanoparticles on the 
I . silica surface, which would lead to an effective surface charge accumulation on the colloids and thus again result in 
electrostatic stabilization. Specifically, the amount of zirconia associated with the microspheres was determined from 
supernatant measurements employing inductively coupled plasma analysis (indicating weak adsorption) as well as 
from scanning angle refiectometry of a silica surface immersed in the nanoparticle solution (indicating no detectable 
adsorption at all). We view the former analysis as more representative, since it involves the actual silica microspheres 
for which the reported phase behavior was determined. 

Since there are only very few approaches for the stabilization of colloidal suspensions, it is of clear importance to 
explore the underlying mechanism of these observations. Not only may a new stabilization technique be of practical 
significance, the binary mixture investigated in Ref. also constitutes an interesting model system for the exploration 
of depletion effects in nonadditive systems 0,|HE1- Lewis and co-workers @] ascribe the initial colloidal stabilization 
to the formation of a "halo" of zirconia particles around the silica colloids, arising from the strong electrostatic repulsion 
between nanoparticles. However, while consistent with zeta-potential measurements that confirm a weak accumulation 
of nanoparticles near the colloidal surface |2( , such halo formation is at variance with the observation of stabilization 
at zirconia volume fractions below 10 -3 , where the average nanoparticle separation greatly exceeds the electrostatic 
screening length A. Indeed, a rather high concentration of nitric acid is required to reach the isoelectric point of silica, 
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resulting in A w 2nm [3, only one third of the nanoparticle diameter. On the other hand, explanation of the reentrant 
gelation representing the upper boundary of the stability window appears more straightforward and is attributed in 
Ref. 0] to the regular depletion attraction 3] induced by the nanoparticles. While this can not be excluded a priori, 
it must be noted that nonadditivity effects can significantly alter the nature of the depletion interaction. 

In Ref. [8j we have studied this system by means of canonical Monte Carlo simulations. While particle-based 
simulations can explicitly account for the fluctuation and correlation effects that are ignored in many analytical 
methods, it is important to note that the extreme size asymmetry a = <T m i cro /er nano = 100 places the current system 
out of reach for conventional Monte Carlo and molecular dynamics simulations. This was resolved by means of a 
novel, highly efficient Monte Carlo scheme that is capable of simultaneously equilibrating interacting species 

of vastly different sizes. In this article, we present a detailed technical account of the calculations reported in Ref. @ 
and make a comparison with results obtained in an integral-equation study of the same system |12L Il3l Il4j . 

Understanding and prediction of phase behavior and stability of a suspension relies on a fundamental knowledge of 
the effective forces between colloids |l5| , which arise from a combination of direct interactions and indirect interactions 
mediated through the solvent and through other solute particles. It is this potential of mean force that we determine 
in the simulations, as a function of nanoparticle concentration. Our findings confirm the experimental observations 
on a quantitative level and, in combination with the corresponding nanoparticle distributions, clarify the physical 
mechanism of the stabilization. Consequently, our computational approach also has predictive capabilities, allowing 
the stability window to be tuned as a function of nanoparticle charge and microsphere-nanoparticle size asymmetry. 



II. COLLOID— NANOPARTICLE INTERACTIONS 



It is the purpose of our calculation to determine the effective interaction V°^ CIO (r) between a pair of colloidal 
microspheres induced by the nanoparticles. Calculation of VJ°? ro (r) amounts to integrating out the degrees of freedom 
of the nanoparticles in the partition function pH Fn| , 

{Snano} 

where H(r; {r}) is the Hamiltonian describing a system containing two microspheres at a separation r and a set of 
nanoparticles with coordinates {r}. V m i CTO (r) represents the direct microsphere pair potential and f3 indicates the 
inverse temperature (k^T) , with fee Boltzmann's constant and T the absolute temperature. The sum runs over all 
nanoparticle configurations s nan o- In order to avoid colloidal many-body effects, the infinite dilution limit must be 
taken for the microspheres. The required direct interactions between nanoparticles and between nanoparticles and 
microspheres are modeled as pairwise potentials V ncmo (r) and V m . n (r), respectively. Thus, the Hamiltonian can be 
written as 

H(r, {r}) = V micro (r) + V wo (\i ij \)+ J2 ^m-n(|iy|) , (2) 

i,j£{nano} i,j'£{m-n} 

where the first sum runs over all pairs of nanoparticles, the second sum runs over all pairs consisting of a micro- 
sphere and a nanoparticle, and = rj — indicates the pair separation. Since K n i cro (r) factorizes in Eq. ft), the 
nanoparticle-induced pair potential is independent of the direct microsphere interaction. 

In our simulation, both colloidal microspheres and nanoparticles are modeled as hard spheres with diameters 
Cmicro and dnano, respectively. We approximate the electrostatic double- layer interactions Vn an o and V m . n via an 
extension of the DLVO theory to nonidentical particles. The Debye-Hiickel approximation is used to linearize 
the Poisson— Boltzmann equation. Subsequent application of the Derjaguin approximation then leads to the Hogg- 
Healy-Fuerstenau (HHF) equation [l^. For a pair of particles at surface-to-surface distance D, under the condition 
of constant surface potential, this equation is given by 

Vhhf(-D) = le Q e r iT ^ (^! 2 + ^ 2 2 ) 



2 (71 + (72 

2*1*9 , 

In 



*1 2 + *2 2 



1 + exp(— kD) 
1 — exp(— kD) 



ln[l -exp(-2/tD)] j , (3) 

where *j and <7j (i = 1, 2) represent the surface potential and diameter, respectively, of a particle of species i, Sq 
denotes the vacuum permittivity, e r = 80 is the dielectric constant of water, and k is the inverse of the Debye length. 

For the interaction between two zirconia nanoparticles, *i — ^2 = ^nano and o~± = o~2 = <7 nano and Eq. @ reduces 
to (Q, Ch. 12) 

Kano(£>) = £o£r7rcr nano * 2 ano ln[l + exp(-KL>)] . (4) 
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Due to the strong screening conditions in Ref. |2|, we can invoke the approximation kD 3> 1 to further simplify this 
expression to 

Kano(-D) = £o£r^ nano *„ ano exp(-fi;.D) . (5) 

While this last approximation overestimates the potential for kD < 1, i.e., for D < kT 1 = 2nm, we note that at 
such small distances even the validity of the original expression Eq. Q becomes less certain. Furthermore, at these 
distances hydration forces become important. As the hydration forces have a range of l-2nm and decay exponentially 
with a very short decay length |20j . we account for the hydration shell of the zirconia particles by replacing <r nano 
with cr nano = Cuano + 1.5nm solely in Eq. 

Since, under the experimental conditions of Ref. Q , the silica microspheres are negligibly charged and thus have 
^micro = 0, we ignore the electrostatic interaction between silica particles. However, Eq. © shows that an electrostatic 
double- layer interaction can arise between the microspheres and the nanoparticles, even though the former are neutral. 
For large size asymmetry a = <7 m j cro /<j nano 3> 1, Eq. © reduces to 

V m - a (D) = ^eoerTTCnano^nano ln[l - exp(-2/tZ>)] , (6) 

which under the approximation kD 3> 1 can be further reduced to 

V m . n (D) = -ie £r7rcrnano*n ano exp(-2K_D) . (7) 

Three important approximations have been made in the derivation of the HHF equation. First, the linearized 
Poisson-Boltzmann theory following from the Debye-Huckel approximation is appropriate only for low surface po- 
tentials. For moderate to high surface potentials, \tj (i = 1, 2) in Eq. has to be replaced b y a n effective surface 
potential Yi(D) = 4(*j/y|) exp(«;D/2) tanh _1 [exp(-K£>/2) tanh(y|/4)] with y? = e^i/k B T [H ]2l| . As a result, the 
HHF equation overestimates the electrostatic interactions for the experimental value \t nano = 70mV by up to 25% 
at large separations. For small values of kD, the effect is considerably less (cf. Fig. 2 in Ref. HJ). Secondly, the 
application of the Derjaguin approximation makes the HHF equation valid only for small kD. However, a strict 
derivation |2lJ without this approximation shows that the HHF formula can be rendered applicable to all kD by re- 
placing (<7io-2)/(<7i + a 2) with {a\a-i)l{a\ + er 2 + 2D). For the short separations and large size ratio investigated here, 
the difference between these factors becomes negligible. Thirdly, a constant surface potential is assumed in Eq. {HI, 
although charge-regulating boundary conditions are more appropriate for the oxide particles employed in Ref. |2j{. 
However, both conditions have been found to lead to comparable results 22], in which an attraction can arise between 
neutral and charged surfaces j^. It is important to emphasize that these mean- field approaches are employed only 
to set up the direct colloid-nanoparticle and nanoparticle-nanoparticle pair potentials to be employed in the Monte 
Carlo simulations. The actual simulations explicitly take into account correlation and fluctuation effects. In addition, 
as is clear from the use of the HHF equation, the solvent and all ions are modeled implicitly, yielding a dielectric 
medium with a specific electrostatic screening length. While the sheer number of solvent molecules and ions precludes 
their explicit incorporation in the calculation, the strong screening also renders this largely unnecessary. Another 
noteworthy point is that the calculations presented here only yield the effective colloidal pair potential. Triplet and 
higher-order interactions, which may become important under the nondilute conditions in the experiments, are not 
addressed. 

The presence of an attractive interaction between microspheres and nanoparticles, as found in Eq. 0, is consis- 
tent with the supernatant measurements reported by Tohver et al. Q. Also their zeta-potential measurements Q 
(which exploit the electrophoretic mobility of the microspheres) indicate a certain degree of microsphere-nanoparticlc 
association, since the increase of the zeta potential with increasing nanoparticle concentration requires a cooperative 
motion of colloids and nanoparticles. Accordingly, it might be tempting to dismiss the haloing phenomenon as a a 
standard case of electrostatic stabilization resulting from an effective charge build-up on the colloids. However, as is 
shown below, the situation is considerably more complicated, since the resulting halo is dynamic in nature and leads 
to effective interactions that are attractive, repulsive, or oscillatory, depending on nanoparticle concentration. 

III. SIMULATION DETAILS 

As already noted in Ref. [lfj, traditional computational approaches experience severe difficulties in the treatment 
of this system. Inherent limitations in most computational methodologies greatly restrict the range of accessible size 
ratios in multi-component mixtures because larger species tend to get trapped by the smaller. This ergodicity problem 
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becomes more severe if small particles accumulate around larger ones, as is the case in the present system |2j,|7J. The 
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TABLE I: Computational parameters for the various nanoparticle concentrations, i^ nmo , investigated. All data apply to a 
nanoparticle-microsphere mixture containing two microspheres at a volume fraction m i cro = 0.002, with a colloid-nanoparticle 
size ratio a — 100. AT nano indicates the number of nanoparticles in the simulation cell. (A^ 1 1 1 " ! f or ) and (A^ano™) indicate 
the average number of microspheres and nanoparticles, respectively, in each cluster. Note that each cluster is started from a 
microsphere. Also shown are the total number of clusters constructed and the required CPU time. 



0nano 


■^V liano 


/ !\rclustcr\ 
V^'micro / 


/ 7V rcl lister \ 
\ nano / 


No. of Clusters 


CPU hours 


0.0007 


700,000 


1.0089 


2503 


200,000,000 


1,000 


0.0010 


1,000,000 


1.0089 


2626 


200,000,000 


1,800 


0.0030 


3,000,000 


1.0092 


5791 


50,000,000 


1,500 


0.0050 


5,000,000 


1.0098 


8737 


40,000,000 


3,100 


0.0070 


7,000,000 


1.0110 


11808 


60,000,000 


4,500 



simulation methods that are capable of explicitly simulating both large and small species are limited to size ratios up 
to a = 10 [24J, except for a few cases of hard-sphere mixtures. Thus, an alternative route that is sometimes taken is 
to only simulate the smaller species around a pair of immobile large particles at a fixed separation [25t I2H |2tJ ■ Since 
the effective force can be evaluated only for one separation at a time, a series of such simulations for different large- 
particle distances must be carried out to construct the entire curve for the effective force. Since the large-particle 
concentration must be kept very small to exclude many-body effects, the total particle number typically becomes 
very large if the small species is present at a fixed concentration and a is large. As a result, this method becomes 
impractically slow even for moderate accuracy. 

We overcome this problem by exploiting a novel, highly efficient Monte Carlo scheme ^| , which permits the 
explicit inclusion of interacting species of vastly different sizes. This algorithm is a variant of the geometric cluster 
algorithm originally introduced by Dress and Krauth [s| . The explicit simulation of both species makes it possible to 
calculate the effective interactions directly from the inversion of the particle pair correlation function g(r) with less 
computational effort than previous techniques. In the infinite dilution limit, the potential of mean force follows from 
Vceif) = —k^T\ng{r). We employ only two microspheres in a periodic cubic box at concentration </> m i C ro = 0.002. 
In principle, the potential of mean force must be calculated from simulations in which the nanoparticles are kept 
at constant chemical potential. For efficiency reasons, we instead use the canonical ensemble, in which the number 
of nanoparticles is constant. Since our simulation cell is very large, fluctuations in the number of nanoparticles are 
expected to be small as the microsphere separation varies, and consequently the difference between the two ensembles 
is expected to be very small |27| . In addition, we have explicitly verified that no changes occur in the effective potential 
as the box size is increased further (i.e., </> m i cr o is lowered further). 

In the actual calculations of the potential of mean force, nano is varied from 0.0007 to 0.01, which encompasses 
the entire range over which stabilization is observed in the experiments. Following the experimental conditions @>@]> 
we choose a microsphere diameter a m i cro = 0.6/^m and a nanoparticle diameter Cnano = 6nm. Furthermore, we set 
$ nano = 70mV and K<7 na no = 3.0, except when explicitly indicated otherwise. Since the colloid-nanoparticle size 
ratio is a = 100, the resulting number of nanoparticles reaches iV nano — 7 x 10 6 for llano = 0.007 (for nano = 0.01, 
where N nano = 10 7 , our calculations are only accurate enough to locate the minimum in the effective potential). The 
corresponding computational effort is summarized in Tabled 

A crucial ingredient to accelerate the simulations is the use of the cell index method with linked lists. Here, two 
independent sets of linked cells are employed, with cell sizes based upon the respective interaction cutoffs used. In 
addition, a direct mapping is used to compute the potential V m . n between species that reside in different cell sets |ll|. 

IV. RESULTS AND DISCUSSION 
A. Nanoparticle adsorption 

In order to investigate the validity of the interactions derived in Sec.|nl we first carry out Monte Carlo simulations 
using the approximate expressions in Eqs. JSJ and Q). The interactions are cut off at r" ano = 5a nano for V nano (r) and 
r™~ n = 60<7 nano for V m . n (r). Unlike in the calculations of the potential of mean force (presented below), here we wish 
to explicitly account for colloidal many-body effects and thus use a system containing 50 microspheres at the same 
volume fraction as the experiments, Amicro = 0.1. Up to 5 x 10 6 nanoparticles are present. 

As illustrated in Fig. the nanoparticle distribution p na no(-D) around a microsphere (expressed as a function of 
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FIG. 1: Log-linear plot of the normalized density profile p n ano(-D)//Onano(oo) of nanoparticles around a microsphere as a function 
of the colloid-nanoparticle surface separation D, for various nanoparticle concentrations </f> nan o- The accumulation near contact 
(D = 0) results from the attractive potential. The minimum that develops with increasing n ano is used to define the extent of 
the halo. 
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FIG. 2: Nanoparticle adsorption per colloidal microsphere as a function of nanoparticle concentration, at fixed microsphere 
concentration, m icro = 0.1. The agreement between experiment (solid squares, Ref. Q) and simulation data (open squares) 
is quite reasonable, in view of the approximate potentials employed. The highest experimental concentration, 3650mg/l, 
corresponds to na no = 10 -3 . 



the microsphere-nanoparticle surface separation D) has a very strong peak near contact, indicating the accumulation 
of nanoparticles around the colloids. The thickness of this accumulation layer, as determined by the minimum of 
Pnano(-D), is approximately <7 nano . For low volume fractions, there is only a negligible difference between the density 
profiles for different na no- However, for higher volume fractions, ^ nano > 0.001, the nanoparticle distribution starts to 
vary, with a more rapid decrease with increasing D and a clear minimum. By integrating p nano {D) from contact to this 
minimum, we compute the amount of adsorbed nanoparticles, shown in Fig.|21along with the experimental findings 2j. 
As can be seen, the agreement with the experiments is fairly good. Both data sets exhibit a linear dependence on n ano 
and a degree of adsorption that is far below 100%. The systematic overestimation in the simulations is compatible 
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^nano ^nano 

FIG. 3: (a) Radial microsphere distribution function g(D) as a function of surface-to-surface separation D, for nanoparticle 
concentrations that correspond to the stable fluid regime, (b) Corresponding effective potential Vcg/k&T between a pair of 
colloidal microspheres. At these low volume fractions, the effective interaction shows a repulsive barrier that prevents colloidal 
aggregation. 

with the high nanoparticle surface potential $ nano = 70mV, which lies somewhat above the regime of validity of 
the Debye-Huckel approximation employed in the HHF equation. At volume fractions above 10 -3 the adsorption 
saturates owing to mutual repulsion between the charged nanoparticles. We conclude that the pair potentials derived 
in Sec. UTI constitute a representative description of the experimental system and capture the essential interactions. 

B. Effective colloidal interactions and colloidal stabilization 

Since the direct interactions have been shown to provide a fairly accurate representation of the experimental system, 
we now proceed to calculate the resulting effective pair potential between the colloids, as a function of nanoparticle 
volume fraction (j) nano - At very low </> n ano = 7 x 1CP 4 the microsphere pair correlation function [Fig. Eta)] shows a 
very strong peak at a surface-to-surface distance D = <7 nan o, corresponding to configurations in which the pair of 
microspheres is separated by a single nanoparticle. However, this peak is followed by a rather broad "exclusion zone," 
corresponding to pair separations that are extremely unlikely to occur. While the short-distance ("bridging") peak 
is only observed for nano = 7 x 1CP 4 , we expect it to be present also for the other concentrations in Fig. |^a), but 
the effective barrier between D = cr nan o and D w 2er nano makes it very unlikely for the system to arrive in such a 
configuration during the course of the simulation. Indeed, the potential of mean force, derived from the microsphere 
pair correlation function and shown in Fig.[3Jb), displays a repulsive barrier of 5k^T at <p nano = 7 x 10~ 4 . If nano is 
increased to 1.0 x 10~ 3 and 3.0 x 10~ 3 , the height of this barrier increases even further. 

Since the van der Waals attraction between microspheres is additive to the nanoparticle-induced interaction (cf. 
Seem, ^ is omitted in the effective potential shown in Fig.^b). While it is strongly attractive at short separations, 
it decays rapidly with increasing microsphere separation and has a strength of approximately only — Ik^T at a surface- 
to-surface separation D = cr nano JtJ. As a result, the repulsive barrier persists in the net pair interaction and is, for 
0nano si 7 x 10 -4 , sufficient to prevent gelation, resulting in kinetic stabilization of the suspension. The corresponding 
threshold agrees remarkably well with the experimentally observed gel-fluid transition near </> nano = 5 x 10~ 4 Q. 

As foreshadowed by the local maximum in g{D) for D « 4er nano at (j) nano — 3.0 x 10~ 3 [Fig. OH a)], the situation 
changes if the nanoparticle volume fraction is raised further. The peak shifts to slightly lower separations and, more 
importantly, increases in height, see Fig. Ufa). This corresponds to the emergence of an effective attraction, which 
reaches a strength of nearly — 3fceT at nano = 10 -2 [Fig. Efb)]. Indeed, since this strength corresponds to the 
typical estimate for the onset of colloidal gelation, we identify nano = 10 -2 as the approximate upper bound of the 
stable fluid region. Again, this is in quite good agreement with the experimentally observed reentrant gelation at 
0nano ~ 5 x 10 -3 , especially in view of the approximate interactions potentials employed. 

A natural explanation for the reentrant behavior might be found in a depletion interaction [j| induced by the 
nanoparticles, as suggested in Ref. 0- However, as shown in the inset of Fig.^fb), the attractive strength increases 
quadratically with </> na no, whereas the classical AO theory for depletion |3j predicts an attractive minimum at con- 
tact {D = 0), with a strength that has a linear concentration dependence. Although the quadratic dependence is 
observed for only less than a decade in 4> nan o, and calculations for larger nanoparticle concentrations are computa- 
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FIG. 4: Counterpart of Fig. 01 but now for higher nanoparticle volume fractions. A peak develops in the radial distribution 
function (a) for surface separations near 4<7 n ano. The corresponding attractive minimum in the effective potential (b) grows 
quadratically with </> na no, as shown in the inset. This attraction is responsible for reentrant gelation of the colloidal microspheres. 



tionally prohibitively expensive, the nonlinearity of the data appears quite pronounced. Since the attractive minimum 
in Fig. Efb) does not occur at contact, one could envisage a situation in which both colloids are decorated with a 
layer of nanoparticles (effectively increasing the colloid diameter) and the remaining nanoparticles induce a deple- 
tion interaction. While this might account for the position of the potential minimum, it does not account for the 
quadratic concentration dependence, which persists even after the corresponding adjustment has been made to the 
nanoparticle concentration employed in the inset of Fig.^b). An alternative explanation for the deviation from the 
linear concentration dependence might be that the strong mutual electrostatic repulsion leads to effective nanoparticle 
concentrations that considerably exceed the bare concentrations, causing the system to leave the regime of validity 
of the AO theory. However, even the effective concentrations are so low that no significant deviations from linearity 
are expected. Thus, we conclude that the effective attractions do not result from a classical depletion mechanism 
alone. Given the nonadditive nature of the interactions, this is not truly surprising. To support this, an examination 
of relevant particle configurations shows no significant nanoparticle depletion between a pair of microspheres at a sep- 
aration D « 4cr nano . It appears therefore plausible that the microsphere attraction arises from correlations between 
the highly charged nanoparticles in the gap between the colloids. 

Since an understanding of the role of size and charge asymmetry between colloids and nanoparticles is of crucial 
importance for potential applications, we have investigated several other parameter combinations. Employing nu- 
merical data from Ref. [23, we have been able to extend the size asymmetry to a = 200 (which involved a system 
with 12 x 10 6 nanoparticles). Figure E^a) displays the effective interaction resulting from a fixed nanoparticle volume 
fraction, nano = 0.003, for colloids of diameter 0.36/im, 0.60//m, and 1.20^m (size ratio a = 60, 100, and 200, re- 
spectively) . For larger a m i cro both the height of the repulsive barrier and the depth of the attractive minimum in the 
effective interaction are enhanced. This is consistent with the experimental observation |2| that both boundaries of 
the stable fluid region were lowered when the size of the colloidal microspheres was increased. The inset of Fig. [SJa) 
shows the same potentials scaled by the size asymmetry. The excellent overlap of all curves shows that the dependence 
on size asymmetry is almost perfectly linear. Thus, the effective interactions follow the Derjaguin approximation, 
according to which the force between two large, identical spheres of diameter a can be expressed in terms of the 
interaction energy of two parallel plates W(D) [2fj|. 

/oo 
( f(z)dz=^W(D), (8) 

where f(z) is the normal force per unit area between two flat surfaces. This relation, which is applicable to any type 
of force law as long as the separation D is much less than the sphere diameter er, indeed predicts that the colloidal 
force (and hence the interaction potential) has a linear dependence on a, as an increase of the colloidal diameter does 
not affect W(D). We note that for actual colloidal stabilization not only the size ratio, but also the absolute size of 
the nanoparticles plays a role, as it determines the range of the repulsive barrier (which must exceed the range of the 
van der Waals interaction). 

Variation of the nanoparticle charge has multiple effects, since it not only alters the nanoparticle repulsion (and 
hence their effective volume fraction and the halo-halo interaction), but also the formation of the halo itself through 
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FIG. 5: (a) Effect of nanoparticle-colloid size asymmetry a on the effective colloidal interactions. At fixed nanoparticle size 
and volume fraction n ano = 0.003, both the strength of the repulsion and the depth of the attractive minimum increase with 
colloid size. As shown in the inset, the effective interactions show a perfectly linear scaling with a, indicating that the Derjaguin 
approximation applies (see text), (b) Effect of nanoparticle charge. All parameters are identical to those in Fig.|3 except for 
the nanoparticle surface potential, which has been lowered from 70mV to 50mV. As a result, much higher nanoparticle volume 
fractions are required to achieve an appreciable colloidal repulsion. 



the colloid-nanoparticle attraction, Eq. J7J. This is also seen by expressing the effective force between a pair colloids 
at surface separation D in terms of the nanoparticle distribution p(r';D), where r' is measured with respect to the 
first colloid, 

F(D) = -Jp(r';D) dVm d ^ r ' l) dr'. (9) 

The nanoparticle charge will affect both p(r'; D) and V m . n { |r'|). Since a systematic variation of the nanoparticle charge 
has not been pursued in experiments, we present simulation results for the zirconia-silica system with a = 100, but 
at a nanoparticle surface potential that has been lowered from $ nano = 70mV to $ nallo = 50mV, see Fig. OJb). 
The effective interactions exhibit the same trend as for higher surface potentials, but at considerably higher volume 
fractions. Indeed, at $ na]10 = 50mV a given effective colloidal repulsion is reached only at a volume fraction that is 
approximately ten times larger than for $ llallo = 70mV. 

Finally, we have performed control simulations in which no attraction is present between the nanoparticles and 
the microspheres. Under these conditions, neither the observed nanoparticle adsorption (Fig. [2Jl nor the effective 
repulsions (Fig. |3J) are recovered, except at nanoparticle concentrations that are much higher than in the experiment 
or if the nanoparticle repulsion has a much longer range than implied by the actual screening length k _1 w cr n ano/3 Q. 



C. Nanoparticle distributions 

Since the effective potentials essentially result from the density distribution of the nanoparticles [cf. Eq. ©], we 
have examined this distribution p(r';D) directly for a pair of colloids at fixed separation D, in order to gain insight 
in the origin of these potentials. Exploiting the symmetry around the z-axis passing through the centers of both 
colloids, we express the nanoparticle positions in cylindrical coordinates (r, (f>, z). A rotational average is performed 
and the r and z coordinates of the nanoparticles are binned. For each bin, the local nanoparticle density is computed. 
The data shown here apply to a smaller size ratio a — 40 (i.e., cr m i cro = 0.24/xm), merely for reasons of clarity. The 
nanoparticle concentration is set to ^ nano = 0.003, resulting in a strongly repulsive effective interaction. Since the 
microspheres are fixed during the simulation, the conventional Metropolis Monte Carlo method is used to simulate 
the movement of the nanoparticles. The averages are obtained over 5,000 independent samples. 

At a surface separation D — 5.0a nano [Fig.|Bfa)] nanoparticles accumulate in the vicinity of the colloidal surface and 
form a uniform, thin halo around each colloid. In this case, the effective interaction is negligibly small. If the colloids 
are forced closer together, the mutual repulsion between the nanoparticles leads to a redistribution of particles in the 
halo. At a separation D = 1.8cr nano [Fig. Efb)] most nanoparticles have been expelled from the gap region between 
the colloids and a depletion zone appears. While reminiscent of the depletion effect in hard-sphere mixtures Q, the 




FIG. 6: [color online] Nanoparticle density distribution around a pair of colloidal particles at a surface-to-surface distance D. (a) 
At larger separations (D = 5.0er nano ), the halos around each colloid are uniform and no appreciable effective interaction occurs, 
(b) Upon closer approach (D — 1.8<T nano ), a depletion of nanoparticles in the gap region is observed, but the effective colloidal 
interaction is strongly repulsive, indicating that energetic terms dominate over entropic terms in the effective interaction. The 
jagged appearance of the distributions results from the binning procedure. 



effective colloidal interaction in this case is strongly repulsive rather than attractive. This indicates that the energy 
penalty resulting from nanoparticle rearrangements dominates over depletion-like entropic effects. As discussed in 
the context of Fig. Efb), at a separation where a weak effective attraction is observed, namely at D w 4.0t7 nano , 
the nanoparticle distribution looks very similar to that depicted in Fig. [Sf a) for D = 5.0t7 nano , i.e., no depletion is 
observed, and detection of significant differences associated with the observed attraction probably requires a spatial 
resolution (bin size) and statistical accuracy that exceeds what is obtained in the current work. 



V. COMPARISON TO INTEGRAL-EQUATION THEORY 

The system studied in this paper has also been investigated by means of integral-equation theory |l2l IT^ | , in which 
the effective colloidal pair potential V^ CTO is computed using the two-com pon ent Ornstein-Zernike (OZ) equations 



supplemented with the hypernetted-chain (HNC) closure. Chavez et al. |12( were the first to demonstrate that, 
under the appropriate conditions, the electrostatic repulsion between nanoparticles alone is sufficient to cause the 
formation of a nanoparticle monolayer on the microsphere surface. Depending on size asymmetry and nanoparticle 
concentration, it is also possible that the nanoparticles are less strongly bound and that a diffuse halo emerges instead 
of a monolayer. Rather than aiming to precisely reproduce the experimental conditions, Ref. |12| focuses on general 
trends. In particular, it is emphasized that, as a general criterion for the formation of a halo or monolayer, the 
microsphere diameter must be compared to the average (bulk) nanoparticle separation. Evidently, this criterion relies 
on a sufficiently strong nanoparticle repulsion. Indeed, all calculations were performed for a specific (large) choice for 
the screening length, KCnano = 0.15, i.e., A « 6.67<7 nano . If seems doubtful that this criterion is applicable for the very 
different conditions of Ref. 0] . 

Subsequently, and coincident with the simulation results 8], another integral-equation study was performed by 
Karanikas and Louis [l^. They find that, for a suitable parameter choice, the effective microsphere interaction 
exhibits a dependence on nanoparticle concentration nan o that is qualitatively compatible with the experimental 
findings. Initially, and just as found already by Chavez et al. 12], the nanoparticles induce a repulsive interaction 
that increases in strength with nanoparticle concentration, resulting in colloidal stabilization at a concentration ^™ c r . 
Upon further increase of 4> nano , an attractive minimum appears in U^f cro , which becomes sufficiently strong to lead 
to reentrant gelation of the colloids at a concentration </>n^o r - 

While the general trends observed in Ref. are certainly encouraging, we note that there are also important 
quantitative discrepancies with the experimental and simulational y| findings. Various calculations in Ref. 
employ a Debye screening length of A — 5.0er nano , much larger than the experimental value A w 0.33er nano and, 
incidentally, rather close to the choice of Chavez et al. [l2| . This increases the mutual repulsion between nanoparticles, 
to the extent where it may alter the mechanism responsible for the observed stabilization. Conversely, in Figure 3 of 
Ref. 13] it is shown that for parameter choices closer to the experimental conditions (e.g., point "A" in this figure, 
corresponding to a screening length A = <r n ano and a nanoparticle contact energy e nano — ft-Ok^T), the lower bound of 
the stability window is shifted to nanoparticle volume fractions that are one to two orders of magnitude larger than 
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FIG. 7: Comparison of the Monte Carlo results of Figs, fflb) and^Ib) (symbols) to HNC integral-equation results (lines) [T^l . 
(a) At low nanoparticle volume fractions, the analytical results reproduce the nanoparticle bridging at D — <r n ano as well as the 
repulsive barrier. Although the height of this barrier is underestimated, it increases with increasing nanoparticle concentration, 
as in the simulations, (b) At higher nanoparticle volume fractions, the integral-equation results reproduce the attractive 
minimum found in Fig.Qfb). The strength of the attraction is approximately twice smaller than in the MC calculations. 



in Refs. UGl, whereas the window is also considerably too narrow, ^no/^nZo ~ 2 instead of ^aSoV^Zo ~ 10 - 
Karanikas and Louis suggest that these differences might stem from the polydisperse nature of the nanoparticles 
employed in the experiments. However, small-angle x-ray scattering measurements |28j have indicated that the 
nanoparticles, while certainly not monodisperse, have a size distribution that is considerably narrower than indicated 
in Ref. Q (2nm < cx nano < 4nm rather than 0.5nm < cx nano < llnm). Nevertheless, these differences are not surprising 
because Karanikas and Louis opted not to include an explicit microsphere-nanoparticle attraction. Accordingly, halo 
formation can only result from nanoparticle repulsion. However, as pointed out in Ref. @ and in Sec. U] of the present 
paper, this mechanism is highly unlikely to be relevant under the conditions of Ref. Q, and can indeed only be made 
to work at either larger screening lengths or higher nanoparticle concentrations. 

While, as shown above and in Ref. y|, the presence of an induced, short-range colloid-nanoparticle attraction can 
resolve this inconsistency, the inclusion of an attractive interaction V m - n yields rather different results for the screening 
length chosen in the integral-equation approach 0] . If V m . n and the nanoparticle repulsion V nano have a comparable 
range, bridging attractions arise that have a detrimental effect on colloidal stability. On the other hand, if V m . n is 
relatively long-ranged compared to Vnano, bridging is avoided, but also the reentrant gelation is no longer observed 
for the entire range of nanoparticle volume fractions investigated |13|. 

It is also of interest to consider instead the predictions of an integral-equation approach for the pair potentials 
employed in the Monte Carlo simulations |g. Figure presents a comparison between the simulation results of Figs. 
I^b) and 01b) and corresponding integral-equation results obtained using the HNC closure |14| . Despite the fact 
that HNC performs best for smoothly varying potentials |29j (whereas the interactions in Sec. [H] vary rapidly at 
short separation), the integral-equation results display semi-quantitative agreement with our simulation results. For 
low (^ n ano [Fig. E^a)], both the bridging attraction at <7 nano and the strong effective repulsion for larger microsphere 
separation are reproduced, although the height of the barrier is systematically underestimated. At </> nan o = 0.0007 
the maximum in V c s is approximately 40% lower than in the MC calculations. For higher nanoparticle concentrations 
[Fig. Efb)] , the secondary minimum near D = 4.0<7 na no is also recovered in the HNC calculations, albeit at a weaker 
strength than the simulations. Overall, we view this as a significant confirmation of the results presented in Ref. jg. 



VI. SUMMARY AND CONCLUSIONS 



We have presented a detailed numerical study of effective interactions between micron-sized silica spheres, induced 
by highly charged zirconia nanoparticles. Our calculations provide an explanation for the nanoparticle haloing phe- 
nomenon discovered by Lewis and co-workers 0, reproducing both colloidal stabilization and reentrant gelation in 
quantitative agreement with the experimental findings. The haloing mechanism for colloidal stabilization is found to 
rely on rather generic features and hence may become of considerable practical importance. Our calculations provide 
explicit guidance for the role of charge and size asymmetry between the colloids and the nanoparticles. The presence 
of weak colloid nanoparticle attractions is found to be a crucial ingredient for stabilization at the very low nanopar- 
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ticle volume fractions and strong screening conditions in the experiments. Integral-equation calculations using the 
hypernetted-chain closure have provided a semi-quantitative confirmation of the effective interactions computed in the 
Monte Carlo simulations. We note the simulations involve not only large size asymmetries, but also extremely large 
numbers of particles. The use of a novel geometric Monte Carlo algorithm |^[l(J,[ll| was therefore indispensable. 
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